2022-09-30 2D Advection-diffusion and waves#

Last time#

  • Implicit Runge-Kutta methods

  • Exploring/discussing tradeoffs

  • SciML benchmarks suite (DifferentialEquations.jl)

Today#

  • FD methods in 2D

  • Cost profile

  • The need for fast algebraic solvers

  • Wave equation and Hamiltonians

  • Symplectic integrators

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end
    
function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))
    x = LinRange(xlims[1], xlims[2], 100)
    y = LinRange(ylims[1], ylims[2], 100)
    heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)
end

struct RKTable
    A::Matrix
    b::Vector
    c::Vector
    function RKTable(A, b)
        s = length(b)
        A = reshape(A, s, s)
        c = vec(sum(A, dims=2))
        new(A, b, c)
    end
end

function rk_stability(z, rk)
    s = length(rk.b)
    1 + z * rk.b' * ((I - z*rk.A) \ ones(s))
end

rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)

function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)
    u = copy(u0)
    t = 0.
    n, s = length(u), length(table.c)
    fY = zeros(n, s)
    thist = [t]
    uhist = [u0]
    while t < tfinal
        tnext = min(t+h, tfinal)
        h = tnext - t
        for i in 1:s
            ti = t + h * table.c[i]
            Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)
            fY[:,i] = f(ti, Yi)
        end
        u += h * fY * table.b
        t = tnext
        push!(thist, t)
        push!(uhist, u)
    end
    thist, hcat(uhist...)
end
ode_rk_explicit (generic function with 1 method)

Extending advection-diffusion to 2D#

1 dimension#

(20)#\[\begin{align} u_t + (- \kappa u_x + w u)_x &= f(x) & \text{ on } \Omega &= (a,b) \\ u(a) &= g_D(a) & u'(b) &= g_N(b) \end{align}\]
  • Cell Peclet number \(\mathrm{Pe}_h = \frac{\lvert w \rvert h}{\kappa}\)

    • \(\mathrm{Pe}_h \lesssim 1\) avoids oscillations

    • \(\mathrm{Pe}_h \gtrsim 1\) is non-stiff for time-dependent model

  • Centered versus upwind for advection

  • Need uniformly bounded \(\kappa \ge \epsilon > 0\)

  • “Strong form” not defined at discontinuities in \(\kappa\)

    • Works okay using divergence form and fluxes at staggered points

2 dimensions#

(21)#\[\begin{align} u_t + \nabla\cdot\big(- \kappa \nabla u + \mathbf{w} u \big) &= f(x,y) & \text{ on } \Omega & \subset \mathbb R^2 \\ u|_{\Gamma_D} &= g_D(x,y) \\ (-\kappa \nabla u + \mathbf w u) \cdot \hat n|_{\Gamma_N} &= g_N(x,y) \end{align}\]
  • \(\Omega\) is some well-connected open set (we will assume simply connected) and the Dirichlet boundary \(\Gamma_D \subset \partial \Omega\) is nonempty.

  • Finite difference methods don’t have an elegant/flexible way to specify boundaries

  • We’ll choose \(\Omega = (-1, 1) \times (-1, 1)\)

On finite difference grids#

  • Non-uniform grids can mesh “special” domains

    • Rare in 3D; overset grids, immersed boundary methods

  • Concept of staggering is complicated/ambiguous

Wesseling 11.4: A boundary-fitted grid around an airfoil.

Time-independent advection-diffusion#

Advection#

\[ \nabla\cdot(\mathbf w u) = \mathbf w \cdot \nabla u + (\nabla\cdot \mathbf w) u\]

If we choose divergence-free \(\mathbf w\), we can use the stencil

\[\begin{split} \mathbf w \cdot \nabla u \approx \frac 1 h \begin{bmatrix} & w_2 & \\ -w_1 & & w_1 \\ & -w_2 & \end{bmatrix} \!:\! \begin{bmatrix} u_{i-1, j+1} & u_{i, j+1} & u_{i+1,j+1} \\ u_{i-1, j} & u_{i, j} & u_{i+1,j} \\ u_{i-1, j-1} & u_{i, j-1} & u_{i+1,j-1} \\ \end{bmatrix} \end{split}\]

Diffusion#

\[ -\nabla\cdot(\kappa \nabla u) = -\kappa \nabla\cdot \nabla u - \nabla\kappa\cdot \nabla u\]
  • When would you trust this decomposition?

  • If we have constant \(\kappa\), we can write

    \[\begin{split} -\kappa \nabla\cdot \nabla u \approx \frac{\kappa}{h^2} \begin{bmatrix} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{bmatrix} \!:\! \begin{bmatrix} u_{i-1, j+1} & u_{i, j+1} & u_{i+1,j+1} \\ u_{i-1, j} & u_{i, j} & u_{i+1,j} \\ u_{i-1, j-1} & u_{i, j-1} & u_{i+1,j-1} \\ \end{bmatrix} \end{split}\]

Advection-diffusion in code#

function advdiff_matrix(n; kappa=1, wind=[1, 1]/sqrt(2))
    h = 2 / n
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    idx(i, j) = (i-1)*n + j
    stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h
    stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2
    for i in 1:n
        for j in 1:n
            if i in [1, n] || j in [1, n]
                push!(rows, idx(i, j))
                push!(cols, idx(i, j))
                push!(vals, 1.)
            else
                append!(rows, repeat([idx(i,j)], 5))
                append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
                append!(vals, stencil_advect + stencil_diffuse)
            end
        end
    end
    sparse(rows, cols, vals)
end
advdiff_matrix (generic function with 1 method)

Spy the matrix#

A = advdiff_matrix(10, wind=[1, 3], kappa=.1)
my_spy(A)
../_images/2022-10-03-2d_15_01.svg
A = advdiff_matrix(20, wind=[0, 0], kappa=.1)
ev = eigvals(Matrix(A))
scatter(real(ev), imag(ev))
../_images/2022-10-03-2d_16_01.svg

Plot a solution#

n = 100
x = LinRange(-1, 1, n)
y = x
f = cos.(pi*x/2) * cos.(pi*y/2)'
heatmap(x, y, f, aspect_ratio=:equal)
../_images/2022-10-03-2d_18_01.svg
A = advdiff_matrix(n, wind=[3,1], kappa=.01)
u = A \ vec(f)
heatmap(x, y, reshape(u, n, n), aspect_ratio=:equal)
../_images/2022-10-03-2d_19_01.svg
  • What happens when advection dominates?

  • As you refine the grid?

Cost breadown and optimization#

using ProfileSVG
function assemble_and_solve(n)
    A = advdiff_matrix(n)
    x = LinRange(-1, 1, n)
    f = cos.(pi*x/2) * cos.(pi*x/2)'
    u = A \ vec(f)
end

@profview assemble_and_solve(200)
Profile results in :-1 #15 in task.jl:484 eventloop in eventloop.jl:8 invokelatest in essentials.jl:726 #invokelatest#2 in essentials.jl:729 execute_request in execute_request.jl:67 softscope_include_string in SoftGlobalScope.jl:65 include_string in loading.jl:1428 eval in boot.jl:368 typeinf_ext_toplevel in typeinfer.jl:996 typeinf_ext_toplevel in typeinfer.jl:1000 typeinf_ext in typeinfer.jl:967 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2340 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2340 abstract_eval_statement in abstractinterpretation.jl:1890 abstract_call in abstractinterpretation.jl:1733 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1600 abstract_apply in abstractinterpretation.jl:1339 abstract_call in abstractinterpretation.jl:1766 abstract_call_known in abstractinterpretation.jl:1696 abstract_call_gf_by_type in abstractinterpretation.jl:162 abstract_call_method_with_const_args in abstractinterpretation.jl:850 typeinf in typeinfer.jl:213 _typeinf in typeinfer.jl:230 typeinf_nocycle in abstractinterpretation.jl:2462 typeinf_local in abstractinterpretation.jl:2366 abstract_eval_statement in abstractinterpretation.jl:1890 _typeinf in typeinfer.jl:257 optimize in optimize.jl:504 run_passes in optimize.jl:539 ssa_inlining_pass! in inlining.jl:85 batch_inline! in inlining.jl:620 ir_inline_item! in inlining.jl:433 just_fixup! in ir.jl:1405 _typeinf in typeinfer.jl:257 optimize in optimize.jl:504 run_passes in optimize.jl:542 sroa_pass! in passes.jl:852 simple_dce! in ir.jl:1422 simple_dce! in ir.jl:1427 iterate in iterators.jl:167 iterate in array.jl:898 assemble_and_solve in In[18]:3 advdiff_matrix in In[2]:1 #advdiff_matrix#1 in In[2]:11 vect in array.jl:126 _array_for in array.jl:679 _array_for in array.jl:676 similar in abstractarray.jl:840 similar in abstractarray.jl:841 Array in boot.jl:468 Array in boot.jl:459 #advdiff_matrix#1 in In[2]:16 repeat in abstractarraymath.jl:356 repeat##kw in abstractarraymath.jl:398 #repeat#1 in abstractarraymath.jl:401 repeat_inner_outer in abstractarraymath.jl:459 repeat_outer in abstractarraymath.jl:479 similar in array.jl:377 Array in boot.jl:459 append! in array.jl:1108 _growend! in array.jl:1011 vect in array.jl:126 _array_for in array.jl:679 _array_for in array.jl:676 similar in abstractarray.jl:840 similar in abstractarray.jl:841 Array in boot.jl:468 Array in boot.jl:459 setindex! in array.jl:966 #advdiff_matrix#1 in In[2]:17 append! in array.jl:1108 _growend! in array.jl:1011 append! in array.jl:1109 copyto! in array.jl:322 _copyto_impl! in array.jl:331 unsafe_copyto! in array.jl:289 vect in array.jl:126 _array_for in array.jl:679 _array_for in array.jl:676 similar in abstractarray.jl:840 similar in abstractarray.jl:841 Array in boot.jl:468 Array in boot.jl:459 #advdiff_matrix#1 in In[2]:18 append! in array.jl:1108 _growend! in array.jl:1011 + in arraymath.jl:16 broadcast_preserving_zero_d in broadcast.jl:849 materialize in broadcast.jl:860 copy in broadcast.jl:885 copyto! in broadcast.jl:913 copyto! in broadcast.jl:957 preprocess in broadcast.jl:940 preprocess_args in broadcast.jl:943 preprocess_args in broadcast.jl:944 preprocess in broadcast.jl:941 broadcast_unalias in broadcast.jl:934 unalias in abstractarray.jl:1427 mightalias in abstractarray.jl:1462 _isdisjoint in abstractarray.jl:1469 != in operators.jl:282 == in promotion.jl:477 #advdiff_matrix#1 in In[2]:22 sparse in sparsematrix.jl:1041 sparse in sparsematrix.jl:1045 sparse in sparsematrix.jl:856 sparse! in sparsematrix.jl:953 setindex! in array.jl:966 sparse! in sparsematrix.jl:1015 setindex! in array.jl:966 assemble_and_solve in In[18]:6 \ in linalg.jl:1564 \ in factorization.jl:105 ldiv! in umfpack.jl:676 ldiv! in umfpack.jl:689 _Aq_ldiv_B! in umfpack.jl:706 _AqldivB_kernel! in umfpack.jl:711 solve! in umfpack.jl:417 umfpack_dl_solve in x86_64-linux-gnu.jl:1768 lu in umfpack.jl:195 #lu#1 in umfpack.jl:0 #lu#1 in umfpack.jl:202 umfpack_numeric! in umfpack.jl:381 #umfpack_numeric!#14 in umfpack.jl:383 umfpack_symbolic! in umfpack.jl:369 umfpack_dl_symbolic in x86_64-linux-gnu.jl:1736 #umfpack_numeric!#14 in umfpack.jl:385 umfpack_dl_numeric in x86_64-linux-gnu.jl:1752
┌ Warning: The depth of this graph is 98, exceeding the `maxdepth` (=50).
│ The deeper frames will be truncated.
└ @ ProfileSVG /home/jed/.julia/packages/ProfileSVG/ecSyU/src/ProfileSVG.jl:275
┌ Warning: The depth of this graph is 98, exceeding the `maxdepth` (=50).
│ The deeper frames will be truncated.
└ @ ProfileSVG /home/jed/.julia/packages/ProfileSVG/ecSyU/src/ProfileSVG.jl:275

What’s left?#

  • Symmetric Dirichlet boundary conditions

  • Symmetric Neumann boundary conditions

  • Verification with method of manufactured solutions

  • Non-uniform grids

  • Upwinding for advection-dominated problems

  • Variable coefficients

  • Time-dependent problems

  • Fast algebraic solvers

Gas equations of state#

There are many ways to describe a gas

Name

variable

units

pressure

\(p\)

force/area

density

\(\rho\)

mass/volume

temperature

\(T\)

Kelvin

(specific) internal energy

\(e\)

energy/mass

entropy

\(s\)

energy/Kelvin

Equation of state#

\[ \rho, e \mapsto p, T \]

Ideal gas#

(22)#\[\begin{align} p &= \rho R T & e &= e(T) \end{align}\]
\[ p = (\gamma - 1) \rho e \]
pressure(rho, T) = rho*T

contour(LinRange(0, 2, 30), LinRange(0, 2, 30), pressure, xlabel="\$\\rho\$", ylabel="\$T\$")
../_images/2022-10-03-2d_27_0.svg

Conservation equations#

Mass#

Let \(\mathbf u\) be the fluid velocity. The mass flux (mass/time) moving through an area \(A\) is

\[ \int_A \rho \mathbf u \cdot \hat{\mathbf n} .\]

If mass is conserved in a volume \(V\) with surface \(A\), then the total mass inside the volume must evolve as

\[ \int_V \rho_t = \left( \int_V \rho \right)_t = - \underbrace{\int_A \rho\mathbf u \cdot \hat{\mathbf n}}_{\int_V \nabla\cdot (\rho\mathbf u)},\]

where we have applied the divergence theorem. Dropping the integrals over arbitrary volumes, we have the evolution equation for conservation of mass.

\[ \rho_t + \nabla\cdot (\rho \mathbf u) = 0 \]

Momentum#

The momentum \(\rho \mathbf u\) has a flux that includes

  • convection \(\rho \mathbf u \otimes \mathbf u\)

    • this is saying that each component of momentum is carried along in the vector velocity field

  • pressure \(p I\)

  • viscous \(-\boldsymbol\tau\)

A similar integral principle leads to the momentum equation

\[ (\rho \mathbf u)_t + \nabla\cdot\big[ \rho \mathbf u \otimes \mathbf u + p I - \boldsymbol \tau \big] = 0 \]

Simplifications#

  • Ignore viscous stress tensor \(\boldsymbol \tau\)

  • Ignore energy equation (not yet written) and assume constant temperature

    • \(p = a^2 \rho\) where \(a\) is the acoustic wave speed

\[\begin{split}\begin{pmatrix} \rho \\ \rho \mathbf u \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \rho \mathbf u \\ \rho \mathbf u \otimes \mathbf u + \rho a^2 I \end{bmatrix} = 0 \end{split}\]

Linearization#

Split each state variable into a mean state and a small fluctuation

  • \(\rho = \bar\rho + \tilde\rho\)

  • \(u = \bar u + \tilde u\)

Let \(\widetilde{\rho u} = (\bar\rho + \tilde\rho) (\bar u + \tilde u) - \bar\rho\bar u \approx \tilde \rho \bar u + \bar\rho \tilde u\), where we have dropped the second order term \(\tilde \rho\tilde u\) because both are assumed small.

We consider background state \(\bar u = 0\) and constant \(\bar\rho(x,y,t) = \bar\rho\). Then

\[\begin{split}\begin{pmatrix} \tilde \rho \\ \bar\rho \mathbf{\tilde u} \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \bar\rho \mathbf{\tilde u} \\ \tilde \rho a^2 I \end{bmatrix} = 0 \end{split}\]

Two forms of acoustic wave equation#

Divide the momentum equation through by background density and dropping the tildes yields the standard form.

\[\begin{split}\begin{pmatrix} \rho \\ \mathbf u \end{pmatrix}_t + \nabla\cdot \begin{bmatrix} \bar\rho \mathbf u \\ \rho \frac{a^2}{\bar\rho} I \end{bmatrix} = 0 .\end{split}\]

Examine second equation

\[ \frac{a^2}{\bar\rho} \nabla\cdot\big[ \rho I \big] = \frac{a^2}{\bar\rho} \nabla \rho \]
and thus $$\begin{pmatrix} \rho \ \mathbf u \end{pmatrix}_t +

(23)#\[\begin{bmatrix} & \bar\rho \nabla\cdot \\ \frac{a^2}{\bar\rho} \nabla & \\ \end{bmatrix}\]
(24)#\[\begin{pmatrix} \rho \\ \mathbf u \end{pmatrix}\]

Let’s differentiate the first equation,

\[ \rho_{tt} + \bar\rho\nabla\cdot(\mathbf u_t) = 0\]
and substitute in the second equation
\[ \rho_{tt} = a^2 \nabla\cdot(\nabla \rho)\]

  • Note: we had to assume these derivatives exist!

We can reduce this to a first order system as

\[\begin{split}\begin{pmatrix} \rho \\ \dot \rho \end{pmatrix}_t + \begin{bmatrix} & -I \\ -a^2 \nabla\cdot\nabla & \end{bmatrix} \begin{pmatrix} \rho \\ \dot\rho \end{pmatrix} = 0\end{split}\]

Question#

  • How is the problem size different?

  • What might we be concerned about in choosing the second formulation?

Laplacian in periodic domain#

function laplacian_matrix(n)
    h = 2 / n
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    wrap(i) = (i + n - 1) % n + 1
    idx(i, j) = (wrap(i)-1)*n + wrap(j)
    stencil_diffuse = [-1, -1, 4, -1, -1] / h^2
    for i in 1:n
        for j in 1:n
            append!(rows, repeat([idx(i,j)], 5))
            append!(cols, [idx(i-1,j), idx(i,j-1), idx(i,j), idx(i+1,j), idx(i,j+1)])
            append!(vals, stencil_diffuse)
        end
    end
    sparse(rows, cols, vals)
end
laplacian_matrix (generic function with 1 method)
L = laplacian_matrix(10)
ev = eigvals(Matrix(L))
scatter(real(ev), imag(ev))
../_images/2022-10-03-2d_40_0.svg

Wave operator#

\[\begin{split}\begin{pmatrix} \rho \\ \dot \rho \end{pmatrix}_t = \begin{bmatrix} & I \\ a^2 \nabla\cdot\nabla & \end{bmatrix} \begin{pmatrix} \rho \\ \dot\rho \end{pmatrix}\end{split}\]
function wave_matrix(n; a=1)
    Z = spzeros(n^2, n^2)
    L = laplacian_matrix(n)
    [Z I; -a^2*L Z]
end
wave_matrix(2)
8×8 SparseMatrixCSC{Float64, Int64} with 16 stored entries:
   ⋅     ⋅     ⋅     ⋅   1.0   ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅   1.0   ⋅    ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅    ⋅   1.0   ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅    ⋅    ⋅   1.0
 -4.0   2.0   2.0    ⋅    ⋅    ⋅    ⋅    ⋅ 
  2.0  -4.0    ⋅    2.0   ⋅    ⋅    ⋅    ⋅ 
  2.0    ⋅   -4.0   2.0   ⋅    ⋅    ⋅    ⋅ 
   ⋅    2.0   2.0  -4.0   ⋅    ⋅    ⋅    ⋅ 
A = wave_matrix(8; a=2) * .1
ev = eigvals(Matrix(A))
plot_stability(z -> rk_stability(z, rk4), "RK4", xlims=(-4, 4), ylims=(-4, 4))
scatter!(real(ev), imag(ev), color=:black)
../_images/2022-10-03-2d_43_0.svg

Question: would forward Euler work?#

Example 2D wave solver with RK4#

n = 20
A = wave_matrix(n)
x = LinRange(-1, 1, n+1)[1:end-1]
y = x
rho0 = vec(exp.(-9*((x .+ 1e-4).^2 .+ y'.^2)))
sol0 = vcat(rho0, zero(rho0))
thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.02)
size(solhist)
(800, 51)
@gif for tstep in 1:length(thist)
    rho = solhist[1:n^2, tstep]
    contour(x, y, reshape(rho, n, n), title="\$ t = $(thist[tstep])\$")
end
┌ Info: Saved animation to 
│   fn = /home/jed/cu/numpde/slides/tmp.gif
└ @ Plots /home/jed/.julia/packages/Plots/lW9ll/src/animation.jl:137

Accuracy, conservation of mass with RK4#

thist, solhist = ode_rk_explicit((t, sol) -> A * sol, sol0, h=.05,
    tfinal=1)

tfinal = thist[end]
M = exp(Matrix(A*tfinal))
sol_exact = M * sol0
sol_final = solhist[:, end]
norm(sol_final - sol_exact)
0.020151117484358525
mass = vec(sum(solhist[1:n^2, :], dims=1))
plot(thist[2:end], mass[2:end] .- mass[1])
../_images/2022-10-03-2d_50_0.svg

Conservation of energy with RK4#

Hamiltonian structure#

We can express the total energy for our system as a sum of kinetic and potential energy.

\[H(\rho, \dot\rho) = \underbrace{\frac 1 2 \int_\Omega (\dot\rho)^2}_{\text{kinetic}} + \underbrace{\frac{a^2}{2} \int_\Omega \lVert \nabla \rho \rVert^2}_{\text{potential}}\]

where we identify \(\rho\) as a generalized position and \(\dot\rho\) as generalized momentum. Hamilton’s equations state that the equations of motion are

\[\begin{split} \begin{pmatrix} \rho \\ \dot\rho \end{pmatrix}_t = \begin{bmatrix} \frac{\partial H}{\partial \dot\rho} \\ -\frac{\partial H}{\partial \rho} \end{bmatrix} = \begin{bmatrix} \dot\rho \\ - a^2 L \rho \end{bmatrix} \end{split}\]

where we have used the weak form to associate \(\int \nabla v \cdot \nabla u = v^T L u\).

function energy(sol, n)
    L = laplacian_matrix(n)
    rho = sol[1:end÷2]
    rhodot = sol[end÷2+1:end]
    kinetic = .5 * norm(rhodot)^2
    potential = .5 * rho' * L * rho
    kinetic + potential
end
ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]
plot(thist, ehist)
../_images/2022-10-03-2d_53_0.svg

Velocity Verlet integrator#

function wave_verlet(n, u0; tfinal=1., h=0.1)
    L = laplacian_matrix(n)
    u = copy(u0)
    t = 0.
    thist = [t]
    uhist = [u0]
    irho = 1:n^2
    irhodot = n^2+1:2*n^2
    accel = -L * u[irho]
    while t < tfinal
        tnext = min(t+h, tfinal)
        h = tnext - t
        u[irho] += h * u[irhodot] + h^2/2 * accel
        accel_next = -L * u[irho]
        u[irhodot] += h/2 * (accel + accel_next)
        accel = accel_next
        t = tnext
        push!(thist, t)
        push!(uhist, copy(u))
    end
    thist, hcat(uhist...)
end
wave_verlet (generic function with 1 method)
thist, solhist = wave_verlet(n, sol0, h=.04)
@gif for tstep in 1:length(thist)
    rho = solhist[1:n^2, tstep]
    contour(x, y, reshape(rho, n, n), title="\$ t = $(thist[tstep])\$")
end
┌ Info: Saved animation to 
│   fn = /home/jed/cu/numpde/slides/tmp.gif
└ @ Plots /home/jed/.julia/packages/Plots/lW9ll/src/animation.jl:137

Accuracy and conservation for Verlet#

thist, solhist = wave_verlet(n, sol0, h=.05, tfinal=50)
tfinal = thist[end]
M = exp(Matrix(A*tfinal))
sol_exact = M * sol0
sol_final = solhist[:, end]
@show norm(sol_final - sol_exact)

mass = vec(sum(solhist[1:n^2, :], dims=1))
plot(thist[2:end], mass[2:end] .- mass[1])
norm(sol_final - sol_exact) = 6.86250099252763
../_images/2022-10-03-2d_58_1.svg
ehist = [energy(solhist[:,i], n) for i in 1:length(thist)]
plot(thist, ehist)
../_images/2022-10-03-2d_59_0.svg

Notes on time integrators#

  • We need stability on the imaginary axis for our discretization (and the physical system)

  • If the model is dissipative (e.g., we didn’t make the zero-viscosity assumption), then we need stability in the left half plane.

  • The split form \(\rho, \rho\mathbf u\) form is usually used with (nonlinear) upwinding, and thus will have dissipation.

Runge-Kutta methods#

  • Easy to use, stability region designed for spatial discretization

  • Energy drift generally present

Verlet/leapfrog/Newmark and symplectic integrators#

  • These preserve the “geometry of the Hamiltonian”

    • energy is not exactly conserved, but it doesn’t drift over time

    • such methods are called “symplectic integrators”

  • May not have stability away from the imaginary axis (for dissipation)

  • Most require a generalized position/momentum split, “canonical variables”